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Abstract 

Complex coherent dynamics is present in a wide variety of neural systems. A typical example is the 
voltage transitions between up and down states observed in cortical areas in the brain. In this work, we 
£N) ' study this phenomenon via a biologically motivated stochastic model of up and down transitions. The 

model is constituted by a simple bistable rate model, where the synaptic current is modulated by short- 
term synaptic processes which introduce stochasticity and temporal correlations. A complete analysis of 
our model, both with mean-field approaches and numerical simulations, shows the appearance of complex 
transitions between high (up) and low (down) neural activity states, driven by the synaptic noise, with 
permanence times in the up state distributed according to a power-law. We show that the experimentally 
observed large fluctuation in up and down permanence times can be explained as the result of sufficiently 
noisy dynamical synapses with sufficiently large recovery times. Static synapses cannot account for this 
behavior, nor can dynamical synapses in the absence of noise. 



u 



x> 



> ■ Author Summary 

vo' 

lO ■ The neural activity observed in most cortical areas often presents a highly complex coherent dynamics. A 

prominent example is the transitions between two well-defined voltage levels measured in different cortical 
regions, that is, the up (high activity) and down (low activity) cortical transitions. Although it is well 
known that the duration of up states is highly irregular (ranging from a few milliseconds to seconds) , the 
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origin of such irregularity is still unclear. In this work we propose that the irregularity in the duration 
of the up states emerges as a consequence of the interplay between synaptic stochasticity and short-term 
plasticity mechanisms. We show, employing analytical treatments and numerical simulations, that such 
interplay induces the appearance of power-law distributions of the permanence times in the up state, 
which may explain the irregularity observed in experiments. On the contrary, such behavior cannot be 



$H ' obtained with static synapses, nor dynamic synapses in absence of noise. 

Introduction 

Neural systems, even in the absence of external stimuli, can exhibit a wide variety of coherent collective 
behaviors, as in vivo and in vitro experiments show [TH3]- O ne °f the most prominent examples is the 
spontaneous transition between two different voltage states, namely up and down states, observed in 
simultaneous individual single neuron recordings as well as in local field measures. Such behavior, which 
is generated within the cortex, may provide a framework for neural computations (3], and could also 
coordinate some sleep rhythms into a coherent rhythmic sequence of recurring cortical and thalamocortical 
activities p2[5J[()]. The phenomenon of up and down transitions has been measured in a number of 
situations, such as in the primary visual cortex of anesthetized animals [H|S], during slow- wave sleep [T][HJ 
[6] , in the somatosensory cortex of awake animals [9] , or in slice preparation under different experimental 
protocols [3[TD1[TT], to name a few. The origin of such structured neuronal activity is still unclear, 
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although several studies have shown that both intrinsic cell properties [T^HH] and the high level of 
recurrency present in actual neural circuits [3l ll5[|TB] may contribute to the generation of up and down 
transitions. In particular, the contribution that reverberations in recurrent neural networks may have 
in the appearance of these transitions could depend strongly on synaptic properties. It is known, for 
instance, that excitatory synapses with a slow dynamics (such as synapses mediated by NMDA receptors) 
may play a relevant role in the generation of persistent activity or up cortical states |17j . On the other 
hand, several modeling studies indicate that activity-dependent synaptic mechanisms, such as short-term 
synaptic depression and facilitation, can induce voltage transitions between up and down neural states 
as well Q1CEBH2!]- 

Many crucial points about the understanding of up and down transitions are, however, still lacking. For 
instance, in vivo experiments in the cat visual cortex show that the permanence times in the depolarized 
or up state present a high variability, and can range from a scale of milliseconds to seconds [7] . A similar 
level of irregularity has also been recently found in in vivo recordings of up-down transitions in the rat 
auditory cortex [2T] , as well as in sleep- wake transitions [15[[2"2"Il2"3] , where power-law distributions in the 
duration of wake states have been measured. Such complexity in the time series of the neuron membrane 
potentials remains far to be explained, and could reflect scale invariance in permanence times, which 
could in turn be a (prcliminar) indicative of criticality. In fact, there are many recent studies that have 
shown criticality in different contexts in the brain [241125] . as well as in neural network models which 
present self-organization and criticality properties [2l)H2"5] , and even it has been reported to occur in 
sleep- wake transitions in in vivo conditions |22[|23j . Although it is worth noting that the irregularity of 
the dynamics of up and down states is not a sufficient condition for criticality, a concrete characterization 
of such irregularity may be a convenient starting point for future works on this topic. 

To study in detail the relevant issue of irregular up and down cortical dynamics, we propose in this 
work a minimal model for up and down transitions in neural media. We consider a simple bistable 
rate model whose stable solutions represent two possible voltage states of the mean membrane potential 
of the network. More precisely, such states correspond, respectively, to high and low levels of activity 
in the network (that is, the up and down cortical states). In addition, we consider that the synaptic 
connections between neurons of the network present short-term synaptic depression (STD) mechanisms, 
which introduce temporal correlations, as well as synaptic stochasticity, in the dynamics of the system 
[29H32J . A complete analysis of this simple mathematical model depicts, both numerically and within 
a theoretical probabilistic approach, the appearance of power-law dependences in the distribution of 
permanence times in the up state. Our results show that the appearance of such scale free distributions 
is due to the complex interplay between several factors including synaptic stochasticity and the temporal 
correlations introduced by STD. The emergence of power-law dependences could, indeed, explain the 
high variability in permanence times in the up state suggested by experiments [7[I21], 

Methods 

Our starting point is a bistable rate model, which mimics the dynamics of the electrical activity of a 
population of interconnected excitatory neurons (although it can be easily extended to other situations) 
with two stable levels of activity. The model has the form [33J 

r„^ = -v{t) + v m S[Jx{t)v{t) -6\+ C(t), (1) 

where v(t) is the mean firing rate of the (homogeneous) neural population, v m is the maximum level of 
activity which can be reached by the population (in absence of noise), J(> 0) is the synaptic coupling 
strength in absence of STD, and 9 is the firing threshold of the neurons in the population. The variable 
£(£) is a Gaussian white noise of zero mean and standard deviation S, which takes into account the inner 
stochasticity of the neural population (caused by other sources of uncontrolled noise in the system). 
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The parameter t v is the population time constant, which may be assumed to be around the duration 
of the synaptic current pulse [341135] . For generality purposes, we set t v = 1, and therefore time and 
frequency are given in units of r v and t~ 1 . respectively. The term S{z) = ^[1 + tanh(z)] represents the 
transduction function, which gives the nonlinear effect that the mean postsynaptic current (coming from 
recurrent connections of the neural population) induces in the network mean firing rate. Employing this 
form for S{z), the up and down stable levels of activity correspond to v ~ v m and v ~ 0, respectively. 

On the other hand, the variable x(t) in equation ^ takes into account the dynamical modification 
of the strength of the synaptic connections during short time scales due to high network activity, and it 
is usually named short-term synaptic plasticity. Based on the model proposed in [291136] for short-term 
depression, and following previous studies concerning the dynamics of neural populations [16j . we assume 
that x{t) evolves according to 

d -^= 1 -^-u X (t H t) + ° m , (2) 

at T r r r 

where r r is the characteristic time scale of the STD mechanism, and u is a parameter related with the 
reliability of the synaptic transmission. According to experimental measurements for these parameters 
in the somatosensory cortex of the rat [36], we set tv = 1000 and u = 0.6 unless specified otherwise! 
The last term on the right hand side of equation ((2J is added to the original model in [36] to include 
some level of stochasticity in this, otherwise, deterministic description of synaptic transmission. The 
inclusion of such term constitutes a simple manner of considering the stochasticity due, for instance, 
to the unreliability of synaptic transmission [31][32], the stochastic properties of receptor-transmitter 
interactions [37], the sparse connectivity of cortical circuits [38l|39], or other sources of noise not yet 
considered (see the Discussion Section for more details). The parameter D controls the strength of this 
fluctuating term, and £(£) is a Gaussian white noise with zero mean and variance one. 

Equations |T]) and @ constitute our minimal model of an excitatory neural network with stochastic 
depressing synapses. The simplifications assumed by such model allows to obtain some analytical deriva- 
tions for the quantities of interest, and concretely for the probability distributions of permanence times 
in the up state, denoted by P(T). Bistable systems in the presence of different sources of noise have been 
theoretically studied in detail in many works [40H44] . Here, however, we have employed a probabilistic 
approach which is very appropriate for the computation of the distribution of permanence times. In the 
following, we will derive an approximate expression for P(T) within this approach. First, we obtain the 
potential function and the conditions in which the dynamics of the system is driven by the variable x. 
After that, we compute the probability distribution of ruin times of x(t) which, as we will see, leads to 
the probability distribution of permanence times in the up state, namely P(T). 

A. The potential function 

In order to compute the potential function of the dynamics (11121) (namely $(f, x)) one can see that, 
for realistic values of tv, the dynamics of x is very slow compared to that of v. We therefore can write 
equation (fTJ) as 

i> = -&,*(!/, a:) + C(<) 

1 u m (3) 

®( v > x ) = t; v ( v - v m) - ttt lo S cosh ( Jxv - e ), 
I Ux 

where we have adiabatically eliminated x from the dynamics of v. The extrema of $ are given by the 

solutions of the equation 

v = -^m [1 + tanh( Jxv - 6>)] = j(y) (4) 



1 Assuming a population time constant of t v = 1 ms, which would approximately correspond to the duration of a fast 
synaptic current pulse mediated by AMPA receptors, we obtain r r = 1000 ms, which is within the physiological range 
measured in 1361 . 
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In the following, we choose 9 = JxqUq, with v$ = \v m an d %o = 1/(1 + ut t vq). With this choice, one can 
easily check from equation ^ that the potential becomes symmetric in v around vq when x ~ Xq. 

Equation (U) may have one or three solutions, depending on the slope of the hyperbolic tangent and 
on the value of 9. In order to obtain three solutions of (@| (that is, the bistable regime) the maximal slope 
of the hyperbolic tangent must be large enough, concretely the condition Jxvq > 1 must be fulfilled. In 
addition, the threshold term must be not too small or too large so that f(v) has three crossing points with 
the straight line v rather than one. This last condition can be written, as a first approach, as f(y\) > v\ 
and fiyi) < V2-, where v\^, are the values where the curvature of the hyperbolic tangent is maximal and 
minimal, respectively. The points V\p. can be easily computed from the third derivative of f(y): 



1 — 3tanh (Jxv — JxqVq) 
cosh (Jxv — JxqVo) 



/'» = v^ 8 ' ,:_~ r , v :"" T ,;,T uy - 

By setting f'"(v) = we obtain 



v x tanh 1 ( v / T/3) 

v x ,2 = ± j^- ■ (6) 

X Jx 

Using now these values for i/ 12 , the conditions f{v\) > v\ and f(y-i) < v 2 can be written as 

1 ,___, _i, r-f-, VqXq /—— 1 , 



-»bVV3+-T-tan^ (vV3) < ^^-^0 < voy/l/S - -ptanh _1 (Vl/3), (7) 

which implies that, in order to have one maxima and two minima in <!>(i/,x), the variable x must be in 
the range x\ < x < X2, where 

_ u Xo + 7 tanh -1 (a/173) _ v xo ~ j tanli -1 (yi73) 

Xi = ■ , , X? = '■ , . (o) 

H,(l + v/I73) ^(l-v/173) 

From equation (jHJ , one can see that the range of x that allows to have three extrema in the potential is 

3 

Ax = X2 — X\ = V3 xq tanh -1 (-\/l/3). (9) 

The condition Ax > implies JxqVq > 1.14 which is, therefore, a sufficient condition to obtain a double 
well potentials for some value of x. Assuming that this condition is satisfied, three different shapes for 
the potential function <b{v,x) can be found, as the figure 1A illustrates. When x < x\ the potential 
function presents only one minimum, located around v ~ 0. Similarly, for X2 < x the potential presents 
also a single minimum, but now located around v ~ v m . Finally, for x\ < x < X2 the potential will take 
a double well shape, with the maximum being located around v ~ vq and the minima located around 
v ~ and v ~ ^ TO , respectively. 

It is worth noting that x\ < xq < X2, with xq being the mean value of x. Due to this, if the range 
Ax is small compared with the fluctuations of x, namely a x , the potential function will spend most of 
the time in the regimes x < x\ and X2 < x, with the double well regime appearing only when the system 
tries to jump from one of these regimes to the other (that is, when x ~ xo). A direct consequence of this 
is that the mean firing rate will be basically switching between the up and down states (that is, v ~ 
and v ~ I'm), and that this switching will be driven by the dynamics of x, as the figure 2 illustrates. 
Therefore, one expects that the distribution of permanence times of v in the up (down) state, becomes 
approximately equal to the distribution of permanence times of x in the x > xq (x < Xq) regime, as long 



2 Onc can find, however, a small discrepancy between our approximate prediction and the actual properties of &(i/,x). 
The discrepancy appears because we have assumed that a sufficient condition for the existence of the three fixed point 
solutions of equation (|4jl is that f{u\) > V\ and /(j^) < ^2, and such assumption is only approximately correct. Plotting 
directly the potential as a function of v reveals that the condition to obtain a double well potential for x ~ xo is Jxqvo > 1, 
rather than Jxqvq > 1.14. 
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as Ax <C <J X is satisficdjj. Due to this equivalence, in order to compute P(T) we only need to compute 
the distribution of permanence times of the variable x in the x > xq regime, denoted as P X (T). 

B. Distribution of permanence times 

In order to compute the distribution of permanence times of x in the x > x$ (or x < xq) regime, one 
can assume that the firing rate takes its mean value v ~ vq in equation ([2]). This is a reasonable approach 
since x is much slower than v for realistic values of the parameters. Considering this approach, and after 
the rescaling z = (1 + uT r v$)x — 1, equation ^ can be written as 

^ = -^ + ^(0 do) 

at t t 

which is the equation of the Ornstein-Uhlenbeck (OU) process (see [55] for details), with r = T r /(l+uT r vo) 
being the correlation time and Zq = z(xo) = 0. Therefore, computing the distribution of permanence 
times in the up state for our system is equivalent to obtain the distribution of the so called ruin timesQ 
for the OU process [¥H1HZ]- The strategy employed here to calculate the distribution of ruin times is 
based on the relation between the ruin time and the first passage time, which is the typical time that 
a stochastic process needs to arrive at a certain threshold value when starting from a certain initial 
condition [47j . Because of the symmetry of the OU process, the distribution of ruin times are equivalent 
when considering excursions of the variable z in the z < region or in the z > region. If we consider 
excursions in the z < region, we can set a small positive threshold e near zero (that is, < e <C 1), 
in such a way that the typical ruin time will be approximately equal to the corresponding first passage 
time, as the figure IB illustrates. The excursions in the region z > typically lead to very short first 
passage times (since e is too small) which we will not take into account in our calculations by considering 
only large enough ruin times. 

The first passage time for the OU process with a small threshold e can be performed by using the 
relation 

7>(e,T\0,0)= f dtV(e,T\e,t)p(t), (11) 

Jo 

where V(a,t a \b,tb) is the conditional probability distribution of the OU process, and p(t) is the first 
passage time distribution. This equation can be solved by taking into account the following property of 
the Laplace transformation 

hit) = f dt' f 2 {t - hit') =» his) = his) his), (12) 

Jo 

where /i(s) is the Laplace transform of hit). By solving the Fokkcr-Planck equation associated with 
equation (fi"0|) . one can obtain the conditional probability for the OU process 

1 ' ' y/2nal[l - cxp(-^t/r)] P \ ^ x [l - exp(-2At/r)] J l ' 



3 It should be noted that, since lis a fraction of available neurotransmitters, its value should be kept within the range 
[0, 1]. In practice, this means that the value of <r x must not be too large, so in order to make Ax <g a x one has to restrict 
to Ax small. In the results presented here, x remain in its realistic range of values, and imposing ad hoc restrictions in such 
a way that x is always within the range [0, 1] docs not affect the results obtained here. 

4 If we consider a stochastic process y(t) starting at t = to from y = yo, the ruin time is denned as the interval t± — to, 
where t\ is the time at which y(t) returns to yo f° r the first time. Since y(t) is a stochastic process, the ruin times are 
stochastic quantities which follow a certain probability distribution. 
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where At = ti — 1\ > 0, and o x = D/y/2r being the standard deviation of x. From expression (fT5)l . and 
assuming that r is large cnough_|, one arrives at 

p(£ ' T|0 ' 0) "73^ cxp (-*) 

(14) 
V(e, Tie, t 1 ) ~ 1 exp (_ ^~*') N ) . 

We denote f x {T) = 7>(e,T|0,0) and / 2 (T - *') = £>(e,T|e,t')- Employing the Laplace transformation in 
/i(T) and /^C? 1 — f) the following expressions are obtained 



fi(s) = Jj^exp[- v / e 2 rs/al 

h( S ) = T/y/e*+4 S TCTl 

Now, taking into account the property (|12p in equation (|11|) . the expression for p(s) is 




^ (s) = V w^ C"^ 5 " 7 " 1 ) ■ (16) 



Finally, for small e one can approximate e 2 + 4stct 2 ~ 4stct 2 . With this approximation, one can easily 
perform the inverse Laplace transformation to equation (|16[) and obtain the distribution of first passage 
times for the OU process 



^-0J-"^(-^f)- (17) 

A similar expression may be obtained if one considers more classical derivations of the first passage 
time of the OU process (see, for instance, [IB]). In order to obtain the distribution of ruin times of the 
OU process, one has to consider a small (but positive) value of e, which leads to p(T) ~ T~ 3 / 2 . The 
distribution of ruin times of the variable x, namely P X {T), and therefore, the distribution of permanence 
times in the up state, namely P(T), for our system are also given by 

P(T)~T~ 3/2 , (18) 

which corresponds to a power-law probability distribution for T. 

Summarizing, the three following conditions must be fulfilled to obtain a power-law dependence in 
P(T) with exponent -3/2: 

• Large enough values of r r . With this condition, we ensure that the dynamics of x(t) is much slower 
than that of v{t). 

• Large enough values of D. In particular, we must have D ^$> 2rAx, according to the condition 
Ax <C <J X and the definitions <j x = D/y2r and t = t,/(1 + uT r vo). This condition can be achieved 
even with very small values of D, since Aa; can be arbitrarily small (by increasing J, for instance). 

• The condition Jxqisq > 1 must hold to ensure the existence of two well defined (up-down) states. 

All these conditions may be easily achieved (up to some point) with realistic values of the model 
parameters, indicating that power-law distributions of the permanence times in the up state are plausible 
to be found in actual cortical media. 



5 More precisely, we assume that r ^> At, which is a valid hypothesis since most of the permanence times in the up state 
are much lower than t. 
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Results 

As we have stated in the previous sections, equations ([T][5]) govern the dynamics of our simplified neural 
system. A typical time series of the dynamics of this model, for the case of deterministic synapses (that 
is, D = 0), is depicted in figure 2 A. In this case, the mean firing rate of the population is characterized 
by a periodic switching between up and down states. This type of periodic behavior was already found 
and analyzed in previous theoretical studies [T4lll6l[T8] and yields bimodal histograms for the mean firing 
rate of the neural population (see figure 2B), as the experiments indicate |3]. However, these approaches 
ignore the stochastic nature of synaptic transmission, and other forms of stochasticity at the synaptic 
level, which seem to be crucial for information processing in neural systems [31 , 32, 49 . Considering 
a certain level of synaptic stochasticity in addition to STD in our model, one obtains a qualitatively 
different emergent behavior, as is shown in figure 2C for D = 20. The mean firing rate presents then 
a complex switching between up and down states, and in particular involves a high variability in the 
permanence times in the up state. 

When deterministic synapses are considered (that is, D = 0) the dynamics of the mean firing rate 
becomes quasi periodic, as it was reported in [T5J[THJ[IH] , for instance. This type of dynamics naturally 
leads to exponential distributions for the permanence timeC|. When D is increased, however, the stochas- 
ticity of the synapses leads to the appearance of power-law distributions for the permanence times in 
the up state. This behavior is shown in figure 3A, where low values of D corresponds to exponential 
distributions for P(T), while larger values of D give P(T) ~ T~ 3 / 2 as predicted by our theoretical cal- 
culations. Such power-law distributions may explain the high variability of permanence times in the up 
state, which has been observed in a number of in vivo experiments, such as in the cat visual cortex [7] 
and rat auditory cortex [21], to name a few. Interestingly, similar power-law dependences have been 
observed during sleep-wake transitions in vivo when one measures the distribution of permanence times 
in the wake state [22j[23]. On the other hand, exponential-like distributions, obtained for the case of 
having D = 0, are not able to explain this variability of the duration of up states. 

By looking at the data for D = 20 in figure 3A, one can observe the existence of a small deviation 
of the numerical results (blue points) with respect to the theoretically predicted slope (solid line) for 
very large values of T. Such deviation is due to the fact that the separation of timescales between the 
dynamics of i/(t) and x(t) (a neccessary condition to obtain power-law dependences in P(T)) is only 
approximate when considering realistic values of the parameters (and in particular, realistic values for 
r r ). More precisely, the approximation fails when the activity of the system falls in the occasional periods 
of very long permanence times in the up state (that is, for large enough values of T, comparable with 
r r ). In order to study the effect of the separation of timescales between v(t) and x(t), we have computed 
P(T) for different (increasing) values of t t while keeping fixed values for Ax and D/r r (this can be 
done by properly modifying J and D with r r , respectively). As a consequence of this, the only effect of 
increasing r r will be a clearer separation of timescales between v(t) and x(t). The results are shown in 
figure 3B, where one can see that larger values of i> (that is, clearer separation of timescales) lead to 
a displacement of the effective cut-off towards higher values of T, as expected, and a clearer power-law 
distribution emerges. 

It is worth noting that the appearance of an effective cut-off in T for realistic conditions does not 
represent an unrealistic feature of the model, but rather it constitutes a prediction about the effective 
range of permanence times which are expected to occur in actual neural systems. Indeed, for realistic 
values of the parameters, our results predict permanence times in the up state up to ~ 1000 ms, which 
is the maximum permanence time observed in experimental realizations j7] . Larger permanence times in 
the up state (of about 10 seconds, for instance) should be expected to appear only as a consequence of 



6 More precisely, for D = our model is similar (except for the term CM) to the one analyzed in 1161 . which shows 
periodic oscillations of the network mean firing rate. In the case of our model with D = 0, the term £(i) introduces certain 
level of stochasticity which turns these periodic oscillations into quasi-periodic oscillations. This leads to the exponential 
distributions for the permanence times in the up state. 
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input driven mechanisms (such as persistent activity associated with working memory tasks |50ll51j ). and 
not as a consequence of spontaneous transitions between different voltage levels, which are the matter of 
interest in this work. 

For a better characterization of the dynamics of the system, one can use, for instance, other statistical 
magnitudes such as the autocorrelation function C(t') of is, which can be defined as 

C(t') = {v{t + *>(*) - v{t)v{t')) . (19) 

Here, (• • • ) indicates a temporal average. The autocorrelation function is depicted in figure 4A for the 
case of deterministic depressing synapses (D = 0) and stochastic depressing synapses {D = 20). C(t') 
presents, for D = 0, two well located peaks at t' ~ ±200, which indicates a strong periodicity of the 
time series (as can be seen in figure 2A). On the contrary, the inclusion of a certain level of intrinsic 
stochasticity in the dynamics of x introduces more pronounced temporal correlations in the dynamics of 
the system. This fact reflects the existence of long permanence stays in the up state, which occurs with 
more probability for high enough values of D, as we have already discussed. 

The spectral properties of the dynamics can be analyzed as well, via the power spectrum defined as 

F(f)= J C{t')exp{2nift / )dt / . (20) 

As one could expect, the power spectrum of the case D = presents a pronounced peak around a 
certain frequency, which in the particular case presented in the figure 4B is / ~ 5 • 10~ 3 . The power 
spectrum for higher values of D shows however different properties than the case D = 0. For instance, 
the figure 4B (which considers D = 20) indicates an approximated power-law behavior for the power 
spectrum, F(f) ~ f~@ with j3 ~ 1.7. This scale-free dependence can be understood by considering that, 
if P(T) is algebraic with exponent 7, the corresponding power spectrum becomes also algebraic with 
exponent j3, where the equation 7 + f3 — 3 relates both exponents |44j . In our particular case, since 
7 ~ 1.5, one obtains a theoretical prediction of (3 ~ 1.5 for the exponent of the power spectrum. The 
theoretical relation between P(T) and F(f) exposed above, however, is only valid under the so called 
single interval approximation, which implies that the integration variable t in equation (|20[) is smaller than 
the permanence time T (see [H] for details). This condition does not strictly hold for our system (where 
T ranges over several scales), and therefore it may introduce deviations in the theoretically predicted 
value of /3 (which is around /? ~ 1.5) with respect to the value found in simulations (of around j3 ~ 1.7). 

Besides the level of synaptic stochasticity, i.e. D, other parameters of the model could have an 
important effect on the dynamics as well. The parameter S, for instance, controls the level of stochasticity 
of the dynamics of v, and therefore one should expect that increasing its value could strongly influence the 
probability distribution P(T). This is shown in figure 5 A, where an increase of 5 disrupts the appearance 
of power-law dependences, and exponential distributions appear instead. This change in P(T) is due to 
the fact that high levels of the additive noise d make the system to jump more frequently from one state 
to the other, and therefore long stays in the up state (and thus distributions with long power-law tails) 
rarely occur. 

The parameters involving the dynamics of x also affect the probability distributions P(T). The 
parameter u, for instance, is responsible for the modulation of x via the mean firing rate v (see equation 
([2])), and therefore it can influence both the dynamics of x and v. As one may see in figure 5B, when 
u takes low values a bump in P(T) emerges for high T. Such deviation from the power-law dependence 
indicates that long stays in the up state occur more frequently than in the power-law case. Attending 
at equation @, one can see that an increase of the mean firing rate v decreases the variable x via the 
parameter u. Therefore, if u takes lower values the decrement of x will be smaller. As a consequence, the 
stays of x in the Xo <C x regime (see Methods Section) will last longer, and the stays of the system in the 
up state will also last longer, causing the observed deviation from the power-law tendency. It should be 
noted, however, that the values of u which allow the appearance of power-law dependences in P(T) for 
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our model agree with the values of u measured in actual cortical media where up and down transitions 
are observed [55] , 

We have also analyzed in detail the effect that varying r r has on the probability distribution of 
permanence times. Note that, contrarily to the previous study presented above, we have now varied the 
parameter r r while all the other parameters are kept fixed. This implies that the modification of r r will 
now have an effect on the separation of timcscales between v(t) and x(t), but also on the concrete value 
of Ax and on the amplitude of the noisy term of equation (J2j (namely D/r r ). The results arc shown in 
figure 5C, where one can distinguish three different regimes as a function of the particular value of r r . 
For low r r (red region in the figure), the probability distributions show an exponential decay for large 
permanence times. The reason for this decay is that, for low r r , the variable x does not perform long 
excursions in the region x$ <C x (sec Methods Section), and therefore the probability to have large values 
of T decreases and the power law behavior for P(T) is not obtained. As r r is increased, long excursions 
for x begin to occur, and we obtain a power law behavior P(T) ~ T~ 3 ' 2 (green region in the figure). 
Finally, one can appreciate that, for even larger values of r r (blue region in the figure), the probability 
distribution of permanence times in the up state presents a power law dependence P(T) ~ T~ l( - Tr ' with 
7(r r ) > 3/2, being an increasing function of r r . Such dependence can not be explained by our previous 
theoretical predictions, based in the assumption that the system is in the bistable regime, and deserves 
a detailed analysis which will be exposed in the next section. 

Further analysis 

In the Methods Section, we established several conditions which had to be fulfilled in order to obtain power 
law dependences for P(T). In particular, our previous analysis indicates that the condition JxqVq > 1 
must hold in order to have a potential function $(^, x) with three extrema (bistable regime). However, as 
we will see in the following, power law expressions for P(T) may appear even if the potential function has 
only one extrcmum in v (concretely, one minimum), although the origin of such power law distributions 
is different from the one considered in previous sections, as we will see. 

When JxqUq < 1 (which occurs for J <C 1 or r r ^> I, for instance), the potential function &(i>, x) 
has only one minimum in is, whose location strongly depends on x. An approximated expression for the 
location of this minimum as a function of x can be obtained by expanding the hyperbolic tangent of 
the fixed point expression of v(t) (see equation (j4j)) around its argument (which is small in this limit), 
yielding 



v„ 



vq{\ - JxqVq) + {Jvq - J^x i/q)x + 3 A vlx\ (21) 



where v m in is the value of v which corresponds to the minimum of the potential function. Therefore as x 
varies around xq , the location of the minimum of the potential v m in also varies in the same way around 
vq. As an example, time series of both v and x are shown in figure 6 A for a given set of parameters 
which satisfies Jxqisq < I. In this time series, the variable v fluctuates around the value v m im which is 
fully determined by x (that is, the variable v becomes a slave variable of x). The predictions of equation 
(|2I[) agree approximately well with simulations and with the numerical evaluation of the fixed points of 
equation ([1]), as the figure 6B shows. 

Since v behaves now as a stochastic variable which does not present a clear bistable dynamics, the 
numerical computation of the distribution of the permanence times will depend on the exact value of v 
above which the system is considered to be in the up state. As we have seen before, this threshold value 
takes the form r\v m (see caption of figure 3), where usually r\ may take a value between 0.6 and 0.9. While 
the results presented for Jx v > 1 (that is, the bistable regime) are quite robust for different values of 
n, in the regime JxoVo this parameter has indeed some effect on P(T), which indicates the difficulty to 
accurately analyze the up and down dynamics in this case. 

In figure 7A, one observes that the distribution P(T) shows also a power law behavior P(T) ~ T 1 
for r] = 0.75 and different values of D, for a set of parameter values which satisfies JxqVo < 1 (that is the 
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monostable regime) . The concrete value of 7 depends strongly on D and it has also a weaker dependence 
with 77, as the figure 7B illustrates. This type of power-law behavior appearing in the monostable regime 
corresponds to the blue region in figure 5C, as well. 

It is worth noting that actual recordings of up and down transitions does not present a clear distinction 
between up and down states, and several nontrivial methods are commonly employed to discriminate 
between both states [52]. Therefore, the results found for the regime Jx^v^ < 1 could indeed reflect the 
behavior of actual cortical up-down transitions, showing power law dependences in P(T) with 7 > 3/2 
and indicating that the concrete nature of the transitions is a synaptic-driven monostable dynamics. 

For a complete characterization of the model, one can summarize all the observed behaviors in a 
phase plot such as the one presented in figure 8A. A total of four different behaviors can be found in 
the (i>, D) space. The first one concerns the dynamics of v whose permanence times in the up state 
follows an exponential distribution (labeled as "E" in the figure) . If the noise amplitude D is sufficiently 
high, one can increase the value of r r to reach the regime "C", in which the dependence P(T) ~ T~ 1 - 5 
is obtained. By increasing r r even more, the probability distribution P(T) takes the form ~ T -7 , with 
7 > 1.5 (regime denoted by "S"), as we have already seen in figure 6. Finally, we also observe that when 
the depression time scale is not large enough (and D < 3), a regime of quasi-periodic- time series of v 
is obtained, with a well-defined duration of up states (regime denoted by "P"). The lines between the 
different regimes have been obtained by visual inspection of P(T) for different values of r r and D. In 
particular, the regime "P" is characterized by the appearance of a bump in the probability distribution 
for some value of T (which reflects a preferred duration of the up state), and the existence of such bump 
has been used as a criterion to distinguish between regimes "P" and "E". Similarly, we assumed that 
the regimes "C" and "S" correspond to the situation in which a power-law behaviour that extends for 
two decades or more is found for P{T). Such criterion, together with an estimation of the slope of the 
power-law via standard Levenberg-Marquardt fitting algorithms, allows to distinguish between regimes 
"E", "C" and "S". 

It must be clarified, however, that actual up and down cortical transitions present most likely a richer 
repertoire of dynamical regimes than the one obtained with our simplified model. It is known, for instance, 
that attractor neural networks with dynamic synapses may exhibit different dynamics corresponding to 
memory, non-memory and switching regimes |18||19) . In this work, we have extensively explored different 
regimes of switching behavior, and its implications for the up and down dynamics observed in the cortex. 
The memory and non-memory regimes, however, can be also found in our simplified model by assuming 
that D, S — > 0. After taking these limits, the system will be in the memory regime if the potential 
function &(v, x) is bistable, or in the non-memory regime if <3>(v, x) is monostable. 

Discussion 

We have shown that the experimentally observed large fluctuations in up and down permanence times 
can be explained as the result of sufficiently noisy dynamical synapses with sufficiently large recovery 
times. Our study suggests that a power-law distribution for these permanence times may emerge as a 
consequence of these two ingredients. Static synapses cannot account for this behavior, nor can dynamical 
synapses in the absence of noise. 

The origin of up and down cortical transitions is still unclear, although different factors that may influ- 
ence their occurrence have been recently reported. It is known, for instance, that inhibitory GABAcrgic 
currents strongly contribute to the temporal coding and spike timing precision of cortical networks during 
up states of activity [3l[53j[54] . Several modeling studies also show the relevance of inhibitory interneu- 
rons in the generation of many types of oscillations in the brain (see for instance j55j). However, other 
studies indicate that most of the main features of up and down transitions depends strongly on synaptic 
plasticity mechanisms, both of long-term and short-term ones [16l[56], and that the transitions appear 
even in the absence of inhibition |16j . In this work we have made the common assumption that the effects 
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of inhibition can be treated as additive and can be incorporated in the threshold of the neuron. This 
is known to be a valid approximation in mean field neural network analysis, but may fail when precise 
timing and details of the dynamical aspects of the neuron affect the inhibition |57U58j . 

Regarding to synaptic characteristics, recent works show that synaptic fluctuations could have an 
important role in the generation of transitions between up and down states 14,59.60 . Since our model 
introduces stochasticity in the synaptic dynamics in a highly simplified manner, however, the last term 
in equation ([2]) should not be associated only with ureliability in synaptic transmission. Indeed, we have 
assumed that other sources of stochasticity may be contributing to this fluctuating term in the mean-field 
quantities i/(t) and x(t). For instance, it is widely known that connectivity in actual cortical media is 
highly sparse. Such feature implies that, in order to obtain the mean-field quantity x(t), the average over 
synapses must be performed over a number of C synapses, with C ranging over 100 ~ 1000 connections 
per neuron |38| . In this situation, the fluctuations of x(t) would be of order 1/vC, which leads to a range 
of 0.1 ~ 0.03 for the values given above. As we have seen, our results state that a value of D/r r — 0.02 is 
enough to obtain power-law distributions (sec figure 3), which lies within this range. Therefore, topology- 
induced fluctuations constitute an important source of stochasticity which could be responsible of the 
appearance of power-law distributions in P(T). Other sources of stochasticity at synaptic level, such 
as the stochastic properties of receptor-transmitter interactions, may also contribute to the last term of 
equation ([2]). Moreover, the low activity rates typical from cortical media lead to a poor time-averaring 
of the incoming input, and therefore the fluctuations at the postsynaptic level will be large at these 
short-time scales (of the order of the typical synaptic integration time constant). 

On the other hand, the amplitude of the noisy term, D/r r , does not need to be very high to induce the 
appearance of power-law distributions in P{T). As we have stated above, a sparse connectivity already 
induces a level of stochasticity which is within the desired range, for instance. Furthermore, the noisy 
term could even be arbitrarily small: attending to our theoretical predictions, a neccessary condition 
to have power-law distributions is that fluctuations of x(t) must be much larger than Ax (see Methods 
Section). Since Ax may be lowered to arbitrary levels (by increasing J, for instance), even a small noisy 
term in the dynamics of x(t) may induce power-law distributions. 

It is also known that short-term synaptic mechanisms, such as short-term depression and facilitation, 
usually play a role in the efficient processing of information. In particular, they may be relevant in many 
tasks, such as in signal detection and coding [29,61 63 or switching between different activity patterns 
previously stored [T91I64] . However, their role on the transitions between cortical states has been pointed 
out only by a few studies [15lH6l[65j . and their possible effects on the statistics of the transitions, which 
is the focus of our work, have been ignored. To the best of our knowledge, the present study is the first 
one which analyzes, even in a simplified manner, the strong effect of synaptic stochasticity- in a general 
sense- and dynamic synapses in the statistics of the up and down transitions. The possible role of other 
short-term synaptic mechanisms, such as STF, has not been addressed yet and constitutes a interesting 
issue still open. 

In our analysis we assumed that the dynamics is symmetric in the up and down states. This is in 
contradiction with experimental evidences |66j which shows that power-law distributions are obtained for 
permanence times in the up state, while permanence times in the down state are exponentially distributed. 
However, this discrepancy disappears when one considers a more realistic transduction function which 
gives an asymmetric potential for the dynamics, and as a consequence the up-down symmetry is broken. 
More detailed studies considering, for instance, some of the biologically realistic aspects discussed above, 
should be performed to test our predictions. In particular, a more elaborated study considering realistic 
neuron models (such as Hodgkin- Huxley model [67]) and stochastic STD models (see [29]I49| . for instance) 
is necessary, as well as more detailed experimental studies which may confirm our predictions. 

From a general point of view, evidences of criticality have been recently found in an increasing num- 
ber of neural systems, such as in the functional connectivity of the living human brain [24], in critical 
avalanches of neuronal activity [25], or in sleep- wake transitions [23], to name a few. According to the 
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results presented in this work, transitions between up and down cortical states could also present some 
relevant properties typical of systems at criticality. Some of these properties have been already mea- 
sured in experiments, such as a high sensitivity of the system to external stimuli [8], or the presence of 
power-law dependences in the power spectra of the neural dynamics [53) . 

It is worth noting that other kind of probability distributions for P(T), such as a log-normal distri- 
bution, could also satisfactorily explain the irregularity in the up states found in experiments. Our study 
shows the importance of some biophysical factors, such as the neurotransmitter recovery time and the 
inherent synaptic stochasticity, and predicts a power-law dependence on P(T) as a consequence of such 
factors. However, further study is needed to investigate other mechanisms, not taken in account in this 
work, which could influence the permanence times in the up state. In a more general sense, our results 
may proportionate a new perspective of the phenomena of up and down transitions (and a theoretical 
framework) that could serve to conciliate the main experimental findings, and that could help for a deep 
understanding of this complex dynamics of the brain activity. 
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Figure 1. Considerations for the mean-field approach. (A) Potential function $(i/, x), as a 
function of the mean firing rate v and for different values of x. One can appreciate the different regimes 
explained in the main text. Other parameters are J = 1.1 V, r r = 1000, u = 0.6 and v m = 5 • 10~ 3 . (B) 
An Ornstein-Uhlenbeck (OU) process (see equation H0| with r = 1000 and D = 20. A typical return 
event (with return time T) and a first passage event (with first passage time T") are indicated for 
illustrative purposes. For the first passage time, the threshold (depicted as a blue dashed line) was fixed 
to 0.15. 
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Figure 2. Time series showing the dynamics of our system. (A) Time series of the mean firing 
rate of the neural population for deterministic depressing synapses. The temporal evolution of the 
variable x is also plotted for illustration purposes. (B) Histogram of the mean firing rate, which shows 
the existence of two well defined states of activity in v ~ 10~ 3 and v ~ 5 • 10 -3 , corresponding to the 
down and up states respectively. The values of the parameters are 

,/ = 1.2 V, T r = 1000, u = 0.6, D = 0, 5 = 0.3 and v m = 5 • 10~ 3 . (C) Same as (A), but with a certain 
level of intrinsic stochasticity on the dynamics of the synapses (concretely, we set D = 20). The 
two-headed arrow shows a typical interval of permanence in the up state, denoted by T. (D) Same as 
(B), but for D = 20. The other parameters take the same values as in (A) and (B). 
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Figure 3. Probability distributions of permanence times in the up state. (A) Probability 
distribution P(T), obtained with numerical simulations, for different values of the noise strength D. 
One can see that high values of D lead to the appearance of power-law distributions P(T) ~ T~ 1 with 
7 = 3/2, as the mean-field solution predicts. For numerical simulations, we employed time series of 
duration 10 6 and averaged over 100 trials. The values of the other parameters were 
J = 1.1 V, u = 0.6, t,. = 1000, 5 = 0.3 and v m = 5 • 10~ 3 . To compute P(T), we have considered that 
the up state has been reached during a period T (with T > 2) if v > nv m during this period. We set 
n = 0.8. (B) Probability distributions of permanence times in the up state, for different values of r r and 
fixed Ax ~ 0.065 and D/r r = 0.02. In order to fix Ax and D/r r , we have conveniently modified J and 
D, respectively, for each value of r r . We employed time series of duration 10 6 and averaged over 600 
trials. Other parameters are u — 0.04, S = 0.3 and v m = 5 • 10 -3 . 
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Figure 4. Autocorrelation and power spectra. (A) Autocorrelation function of the mean firing 
rate for deterministic (D = 0) and stochastic (D = 20) synapses, in the presence of STD. (B) Power 
spectra of the mean firing rate for the two cases illustrated in (A). For both panels, we have averaged 
over 10 5 time series of duration 10 6 each, and we have fixed J = 1.1 V, u=- 0.6, r,. = 1000, 5 = 0.3 and 
v m = 5 • 10~ 3 . 
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Figure 5. Influence of other parameters of the model. (A) Probability distributions of 

permanence times in the up state, for different values of 6. Other parameters are 

,/ = 1.1 V, u = 0.6, r r = 1000, D = 20 and v m = 5 • 10~ 3 . (B) Same as in (A), but for different values 

of u. The other parameters take the same values as in (A), except for 8 = 0.3. (C) Probability 

distribution P(T) as a function of T and r r . The three different regimes are shown with different colors 

(see main text for details). Other parameters are J = 1.1 V, u = 0.6, D = 20, 5 = 0.3 and 

v m = 5 • 10~ 3 . For all panels, we have averaged over 100 times series of duration 10 6 each. 
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Figure 6. Behavior of the system when the condition Jx$vq < 1 holds. (A) Time series of the 
variables v and x. (B) The same time series, but represented on the x — v plane, illustrates the fact that 
v is a slave variable of x (although some level of inner stochasticity on v is still present). The green line 
corresponds to the approximate expression (f2Tj) , while the blue line is the numerical evaluation of the 
fixed point solutions of v(t) (see equation (J4J). The inset shows the situation in which the system shows 
a bistable dynamics, analyzed in the previous section. (C) The potential function as a function of v for 
different values of x. One can appreciate the existence of only one minimum, whose location is 
controlled by x. (D) Histograms of the mean firing rate of the system for different values of J. For the 
cases showed in this panel, the condition JxqVo < 1 is only satisfied for the case J = 0.55. For all 
panels, u = 0.6, t t = 1000, D = 20, S = 0.3, v m = 5 • 10~ 3 , and J = 0.55 V unless specifically specified. 



Irregular up and down transitions 



22 



A 



B 




?- 



1000 




200 



Figure 7. Statistics of permanence times in the up state for JxqVq < 1. 
distribution of permanence times in the up state in the JxqVq < 1 regime, for n 



(A) Probability 
= 0.75 and different 



values of D. One can see that power law relations P(T) ~ T~ 7 appear. (B) Dependence of 7 with D for 
the conditions presented in (A) . The inset shows the dependence of 7 with the parameter r\ for the case 
D = 200. We have averaged over 100 time series of duration 10 6 each. Other parameters are 
J = 0.05 V, 7v = 1000, u = 0.6, 8 = 0.3 and v m = 5 • 10~ 3 . 
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Figure 8. The different dynamical regimes of the model. (A) Phase plot which shows the 
different behaviors found in our system. These behaviors corresponds to time series of v for which 
permanence times in the up state follow an exponential distribution (E), a power-law distribution 
P{T) ~ T 1 with 7 = 3/2 (C), or a power-law distribution with 7 > 3/2 (S). In addition, a phase with 
a well-defined duration of the up state is found (P). In panel (B) some of these behaviors are depicted. 
From top to bottom one can see situations P, E and C. Other parameters are 
J = 1.1 V, u = 0.6, 5 = 0.3 and v m = 5 ■ 10~ 3 . 



